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Abstract 

The reconstruction of a 3D object or a scene is a classical inverse 
problem in Computer Vision. In the case of a single image this is 
called the Shape-from-Shading (SfS) problem and it is known to be 
ill-posed even in a simplified version like the vertical light source case. 

A huge number of works deals with the orthographic SfS problem based 
on the Lambertian reflectance model, the most common and simplest 
model which leads to an eikonal type equation when the light source is 
on the vertical axis. In this paper we want to study non-Lambertian 
models since they are more realistic and suitable whenever one has 
to deal with different kind of surfaces, rough or specular. We will 
present a unified mathematical formulation of some popular ortho¬ 
graphic non-Lambertian models, considering vertical and oblique light 
directions as well as different viewer positions. These models lead to 
more complex stationary nonlinear partial differential equations of Ha¬ 
milton-Jacobi type which can be regarded as the generalization of the 
classical eikonal equation corresponding to the Lambertian case. How¬ 
ever, all the equations corresponding to the models considered here 
(Oren-Nayar and Phong) have a similar structure so we can look for 
weak solutions to this class in the viscosity solution framework. Via 
this unified approach, we are able to develop a semi-Lagrangian ap¬ 
proximation scheme for the Oren-Nayar and the Phong model and to 
prove a general convergence result. Numerical simulations on synthetic 
and real images will illustrate the effectiveness of this approach and the 
main features of the scheme, also comparing the results with previous 
results in the literature. 
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1 Introduction 


The 3D reconstruction of an object starting from one or more images is a 
very interesting inverse problem with many applications. In fact, this prob¬ 
lem appears in various fields which range from the digitization of curved 
documents m to the reconstruction of archaeological artifacts m- More 
recently, other applications have been considered in astronomy to obtain 
a characterization of properties of planets and other astronomical entities 
[MlEaEI] and in security where the same problem has been applied to the 
facial recognition of individuals. 

In real applications, several light sources can appear in the environment 
and the object surfaces represented in the scene can have different reflection 
properties because they are made by different materials, so it would be hard 
to imagine a scene which can satisfy the classical assumptions of the 3D 
reconstruction models. In particular, the typical Lambertian assumption 
often used in the literature has to be weakened. Moreover, despite the fact 
that the formulation of the Shape-from-Shading problem is rather simple 
for a single light source and under Lambertian assumptions, its solution is 
hard and requires rather technical mathematical tools as the use of weak 
solutions to nonlinear partial differential equations (PDEs). From the nu¬ 
merical point of view the accurate approximation of non regular solutions to 
these nonlinear PDEs is still a challenging problem. In this paper we want 
to make a step forward in the direction of a mathematical formulation of 
non-Lambertian models in the case of orthographic projection with a sin¬ 
gle light source located far from the surface. In this simplified framework, 
we present a unified approach to two popular models for non-Lambertian 
surfaces proposed by Oren-Nayar [IBl SH SZl |50] and by Phong [5l]. We 
will consider light sources placed in oblique directions with respect to the 
surface and we will use that unified formulation to develop a general numeri¬ 
cal approximation scheme which is able to solve the corresponding nonlinear 
partial differential equations arising in the mathematical description of these 
models. 

To better understand the contribution of this paper, let us start from the 
classical SfS problem where the goal is to reconstruct the surface from a sin¬ 
gle image. In mathematical terms, given the shading informations contained 
in a single two-dimensional gray level digital image /(x), where x := 
we look for a surface 2 ; = i^(x) that corresponds to its shape (hence the 
name Shape from Shading). This problem is described in general by the 
image irradiance equation introduced by Bruss [8] 

/(x) = i?(N(x)), (1) 

where the normalized brightness of the given grey-value image /(x) is put 
in relation with the function i?(N(x)) that represents the reflectance map 
giving the value of the light reflection on the surface as a function of its 
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orientation (i.e., of the normal N(x)) at each point (x,?i(x)). Depending on 
how we describe the function i?, different reffection models are determined. 
In the literature, the most common representation of R is based on the 
Lambertian model (the L-model in the sequel) which takes into account 
only the angle between the outgoing normal to the surface N(x) and the 
light source u;, that is 

I(x) = 7 £,(x)N(x)-CJ, (2) 

where • denotes the standard scalar product between vectors and 7 d(x) in¬ 
dicates the diffuse albedo, i.e. the diffuse reffectivity or reffecting power of a 
surface. It is the ratio of reffected radiation from the surface to incident ra¬ 
diance upon it. Its dimensionless nature is expressed as a percentage and it 
is measured on a scale from zero for no reffection of a perfectly black surface 
to 1 for perfect reffection of a white surface. The data are the grey-value 
image /(x), the direction of the light source represented by the unit vector 
(jj and the albedo 7 d(x). The light source a; is a unit vector, hence |a;| = 1. 
In the simple case of a vertical light source, that is when the light source is 
in the direction of the vertical axis, this gives rise to an eikonal equation. 
Several questions arise, even in the simple case: is a single image sufficient 
to determine the surface? If not, which set of additional informations is nec¬ 
essary to have uniqueness? How can we compute an approximate solution? 
Is the approximation accurate? It is well known that for Lambertian sur¬ 
faces there is no uniqueness and other informations are necessary to select 
a unique surface (e.g. the height at each point of local maximum for /(x)). 
However, rather accurate schemes for the classical eikonal equation are now 
available for the approximation. Despite its simplicity, the Lambertian as¬ 
sumption is very strong and does not match with many real situations that 
is why we consider in this paper some non-Lambertian models trying to give 
a unified mathematical formulation for these models. 

In order to set this paper into a mathematical perspective, we should men¬ 
tion that the pioneering work of Horn [281 EH] and his activity with his 
collaborators at MIT [301 EH produced the first formulation of the Shape 
from Shading (SfS) problem via a partial differential equation (PDE) and 
a variational problem. These works have inspired many other contributions 
in this research area as one can see looking at the extensive list of refer¬ 
ences in the two surveys jESllIO]. Several approaches have been proposed, 
we can group them in two big classes (see the surveys [ESI Hi): methods 
based on the resolution of PDEs and optimization methods based on a vari¬ 
ational approximation. In the first group the unknown is directly the height 
of the surface z = i^(x), one can find here rather old papers based on the 
method of characteristics miEniEaisiiisiEiEHi where one typically looks 
for classical solutions. More recently, other contributions were developed 
in the framework of weak solutions in the viscosity sense starting from the 
seminal paper by Rouy and Tourin m and, one year later, by Lions-Rouy 
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and Tourin |1D] (see e.g. [33l El El |22l [101 EH SI ISSl (Ml 1231 [53]). The 
second group contains the contribution based on minimization methods for 
the variational problem where the unknown are the partial derivatives of 
the surface, p — Ux and q — Uy (the so-called normal vector field. See 
e.g. [301 [ini EZl [Ml [Ml [321 El [fe]). It is important to note that in this 
approach one has to couple the minimization step to compute the normal 
field with a local reconstruction for u which is based usually on a path in¬ 
tegration. This necessary step has also been addressed by several authors 
(see m and references therein). We should also mention that a continuous 
effort has been made by the scientific community to take into account more 
realistic reflectance models [3 E Sa Ea EH, different scenarios including 
perspective camera projection mi EH El E3 EH and/or multiple im¬ 
ages of the same object [83l |8l|. The images can be taken from the same 
point of view but with different light sources as in the photometric stereo 
method ISD Ezi ia sa US] or from different points of view but with the 
same light source as in stereo vision m- Recent works have considered 
more complicated scenarios, e.g. the case when light source is not at the 
optical center under perspective camera projection [35|- R is possible to 
consider in addition other supplementary issues, as the estimation of the 
albedo [sa la Es iM] or of the direction of the light source that are usu¬ 
ally considered known quantities for the model but in practice are hardly 
available for real images. The role of boundary conditions which have to be 
coupled with the PDE is also a hard task. Depending on what we know, the 
model has to be adapted leading to a calibrated or uncalibrated problem 
(see miiHsiEaEzi for more details). In this work we will assume that the 
albedo and the light source direction are given. 

Regarding the modeling of non-Lambertian surfaces we also want to men¬ 
tion the important contribution of Ahmed and Farag [2]. These authors have 
adopted the modeling for SfS proposed by Prados and Faugeras [S3IS2] us¬ 
ing a perspective projection where the light source is assumed to be located 
at the optical center of the camera instead at infinity and the light illumi¬ 
nation is attenuated by a term 1/r^ (r represents here the distance between 
the light source and the surface). They have derived the Hamilton-Jacobi 
(HJ) equation corresponding to the Oren-Nayar model, developed an ap¬ 
proximation via the Lax-Friedrichs sweeping method. They gave there an 
experimental evidence that the non-Lambertian model seems to resolve the 
classical concave/convex ambiguity in the perspective case if one includes 
the attenuation term 1/r^. In [I] they extended their approach for various 
image conditions under orthographic and perspective projection, comparing 
their results for the orthographic L-model shown in [63] and in [85]. Finally, 
we also want to mention the paper by Ragheb and Hancock [5H] where they 
treat a non-Lambertian model via a variational approach, investigating the 
reflectance models described by Wolff and by Oren and Nayar iHniEsiisnj. 
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Our Contribution. 

In this paper we will adopt the PDE approach in the framework of weak 
solutions for Hamilton-Jacobi equations. As we said, we will focus our at¬ 
tention on a couple of non-Lambertian reflectance models: the Oren-Nayar 
and the Phong models liHiiiniiiTiEniEi]. Both models are considered for 
an orthographic projection and a single light source at inflnity, so no atten¬ 
uation term is considered here. We are able to write the Hamilton-Jacobi 
equations in the same flxed point form useful for their analysis and approxi¬ 
mation and, using the exponential change of variable introduced by Kruzkov 
in [39], we obtain natural upper bound for the solution. Moreover, we pro¬ 
pose a semi-Lagrangian approximation scheme which can be applied to both 
the models, prove a convergence result for our scheme that can be applied 
to this class of Hamilton-Jacobi equations, hence to both non-Lambertian 
models. Numerical comparisons will show that our approach is more accu¬ 
rate also for the 3D reconstructions of non-smooth surfaces. 

A similar formulation for the Lambertian SfS problem with oblique light 
direction has been studied in [23] and here is extended to non-Lambertian 
models. We have reported some preliminary results just for the Oren-Nayar 
model in m 

Organization of the paper. 

The paper is organized as follows. After a formulation of the general model 
presented in Section we present the SfS models starting from the classi¬ 
cal Lambertian model (Section |^. In Sections and we will give details 
on the construction of the nonlinear partial differential equation which cor¬ 
responds respectively to the Oren-Nayar and the Phong models. Despite 
the differences appearing in these non-Lambertian models, we will be able 
to present them in a unifled framework showing that the Hamilton-Jacobi 
equations for all the above models share a common structure. Moreover, the 
Hamiltonian appearing in these equations will always be convex in the gra¬ 
dient Vu. Then, in Section]^ we will introduce our general approximation 
scheme which can be applied to solve this class of problems. In Sectionwe 
will apply our approximation to a series of benchmarks based on synthetic 
and real images. We will discuss some issue like accuracy, efficiency and the 
capability to obtain the maximal solution showing that the semi-Lagrangian 
approximation is rather effective even for real images where several param¬ 
eters are unknown. Finally, in the last section we will give a summary of 
the contributions of this work with some final comments and future research 
directions. 

2 Formulation of the general model 

We fix a camera in a three-dimensional coordinate system {Oxyz) in such a 
way that Oxy coincides with the image plane and Oz with the optical axis. 
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Let a; = (001^002^003) = (00^003) G (with 003 > 0) be the unit vector that 
represents the direction of the light source (the vector points from the object 
surface to the light source); let /(x) be the function that measures the gray- 
level of the input image at the point x := (x^y). /(x) is the datum in the 
model since it is measured at each pixel of the image, for example in terms 
of a greylevel (from 0 to 255). In order to construct a continuous model, 
we will assume that /(x) takes real values in the interval [0,1], defined in a 
compact domain Q called “reconstruction domain” (with C open set), 
/ : 12 ^ [0,1], where the points with a value of 0 are the dark point (blacks), 
while those with a value of 1 correspond to a completely reflection of the 
light (white dots, with a maximum reflection). 

We consider the following assumptions: 

Al. there is a single light source placed at infinity in the direction 00 (the 
light rays are, therefore, parallel to each other); 

A2. the observer’s eye is placed at an infinite distance from the object you 
are looking at (i.e. there is no perspective deformation); 

A3, there are no autoreflections on the surface. 


In addition to these assumptions, there are other hypothesis that depend on 
the different reflectance models (we will see them in the description of the 
individual models). 

Being valid the assumption (A2) of orthographic projection, the visible part 
of the scene is a graph 2 : = u(x) and the unit normal to the regular surface 
at the point corresponding to x is given by: 


N(x) 


n(x) 

n(x)| 


(-V^(x),l) 
yi + |Vn(x)|2’ 


(3) 


where n(x) is the outgoing normal vector. 

We assume that the surface is standing on a flat background so the height 
function, which is the unknown of the problem, will be non negative, u : 
fl [0,oc). We will denote by fl the region inside the silhouette and we 
will assume (just for technical reasons) that fl is an open and bounded sub¬ 
set of M? (see Fig. [^. It is well known that the SfS problem is described 
by the image irradiance equation 0 and depending on how we describe the 
function R different reflection models are determined. We describe below 
some of them. To this end, it would be useful to introduce a representation 
of the brightness function / (x) where we can distinguish different terms rep¬ 
resenting the contribution of ambient, diffuse reflected and specular reflected 
light. We will write then 


I(x) = /ca/^(x) + knlDi^) + (4) 
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Figure 1: An object on a flat background: Q indicates the region inside the 
silhouette, dQ the boundary of it. 


where Ia{^)-> Id{^) and /s'(x) are respectively the above mentioned compo¬ 
nents and kA^ ko and ks indicate the percentages of these components such 
that their sum is equal to 1 (we do not consider absorption phenomena). 
Note that the diffuse or specular albedo is inside the definition of /^(x) or 
/s'(x), respectively. In the sequel, we will always consider /(x) normalized in 
[0,1]. This will allow to switch on and off the different contributions depend¬ 
ing on the model. Let us note that the ambient light term /a(x) represents 
light everywhere in a given scene. In the whole paper we will consider it 
as a constant and we will neglect its contribution fixing kA — 0. Moreover, 
for all the models presented below we will suppose uniform diffuse and/or 
specular albedo and we will put them equal to 1, that is all the points of 
the surface reflect completely the light that hits them. We will omit them 
in what follows. As we will see in the following sections, the intensity of 
diffusely reflected light in each direction is proportional to the cosine of the 
angle 9i between surface normal and light source direction, without taking 
into account the point of view of the observer, but another diffuse model 
(the Oren-Nayar model) will consider it in addition. The amount of specu¬ 
lar reflected light towards the viewer is proportional to (cos^^)"^, where 0s 
is the angle between the ideal (mirror) reflection direction of the incoming 
light and the viewer direction, a being a constant modelling the specularity 
of the material. In this way we have a more general model and, dropping 
the ambient and specular component, we retrieve the Lambertian reflection 
as a special case. 

3 The Lambertian model (L—model) 

For a Lambertian surface^ which generates a purely diffuse model, the specu¬ 
lar component does not exist, then in @ we have just the diffuse component 
Id on the right side. Lambertian shading is view independent, hence the 
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(5) 


irradiance equation 0 becomes 

/(x) = N(x) • a?. 

Under these assumptions, the orthographic SfS problem consists in deter¬ 
mining the function ^ M that satisfies the equation 0 . The unit 

vector uj and the function I (x) are the only quantities known. 

For Lambertian surfaces [201EU , just considering an orthographic projection 
of the scene, it is possible to model the SfS problem via a nonlinear PDF 
of the first order which describes the relation between the surface i/(x) (our 
unknown) and the brightness function /(x). In fact, recalling the definition 
of the unit normal to a graph given in Q, we can write ^ as 

/(x) ^1 + |Vi/(x)p + u; • Vi^(x) — 6^3 = 0, in ( 6 ) 

where uj = This is an Hamilton-Jacobi type equation which does 

not admit in general a regular solution. It is known that the mathematical 
framework to describe its weak solutions is the theory of viscosity solutions 
as in |i 0 ] . 

The vertical light case. 

If we choose a; = (0,0,1), the equation Q becomes the so-called “eikonal 
equation”: 

|Vi/(x)| = /(x) for X G iJ, (7) 

where 

The points x G where /(x) assumes maximum value correspond to the 
case in which a; and N(x) have the same direction: these points are usually 
called “singular points”. 

In order to make the problem well-posed, we need to add boundary condi¬ 
tions to the equations 0 or 0 : they can require the value of the solution 
u (Dirichlet boundary conditions type), or the value of its normal derivative 
(Neumann boundary conditions), or an equation that must be satisfied on 
the boundary (the so-called boundary conditions “state constraint”). In this 
paper, we consider Dirichlet boundary conditions equal to zero assuming a 
surface on a flat background 

u{-k) = 0, for X G 9D, (9) 

but a second possibility of the same type occurs when it is known the value 
of u on the boundary, which leads to the more general condition 


?x(x)=^(x), for X G 9D. 


(10) 





Unfortunately, adding a boundary condition to the PDE that describes the 
SfS model is not enough to obtain a unique solution because of the con¬ 
cave/convex ambiguity. In fact, the Dirichlet problem Q-(IO) can have 
several weak solutions in the viscosity sense and also several classical so¬ 
lutions due to this ambiguity (see m)- As an example, all the surfaces 
represented in Fig. are viscosity solutions of the same problem 0 -© 
which is a particular case of 0 -([Tol) (in fact the equation is \u'\ = —2x 
with homogenous Dirichlet boundary condition). The solution represented 
in Fig. [^a is the maximal solution and is smooth. All the non-smooth a.e. 
solutions, which can be obtained by a reflection with respect to a horizontal 
axis, are still admissible weak solutions (see Fig. [^b). In this example, 
the lack of uniqueness of the viscosity solution is due to the existence of 
a singular point where the right hand side of Q vanishes. An additional 
effort is then needed to deflne which is the preferable solution since the lack 
of uniqueness is also a big drawback when trying to compute a numerical 
solution. In order to circumvent these difficulties, the problem is usually 
solved by adding some information such as height at each singular point 

m- 




Figure 2: Illustration of the concave/convex ambiguity: (a) maximal solu¬ 
tion and (b) a.e. solutions giving the same image. Figure adapted from 

m- 

For analytical and numerical reasons it is useful to introduce the expo¬ 
nential Kruzkov change of variable [39] jxv{'x.) = 1 — In fact, setting 

the problem in the new variable v we will have values in [0, l//i] instead of 
[0, oc) as the original variable u so an upper bound will be easy to And. Note 
that /i is a free positive parameter which does not have a speciflc physical 
meaning in the SfS problem. However, it can play an important role also 
in our convergence proof as we will see later (see the remark following the 
end of Theorem |6.1[ ). Assuming that the surface is standing on a flat back¬ 
ground and following |23| , we can write 0 and 0 in a fixed point form in 
the new variable v. To this end let us define x 9^3(0,1) ^ and 

y : X dBziO, 1) X [0,1] R as 

b^(x, a) := — (I(x)ai -6ai,I(x)a2 -^2), (11) 
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/■^(x, a, i;(x)) := - /xt>(x)) + 1 

CJS 


(12) 


and let denote the unit ball in We obtain 


Lambertian Model 



(x, 'y(x), V'L’) for X G ri, 
for X G dLt 


(13) 


where 


T^(-) := min {b^(x, a) • V'?;(x) + /^(x, a, '?;(x))}. 


It is important to note for the sequel that the structure of the above first 
order Hamilton-Jacobi equation is similar to that related to the dynamic 
programming approach in control theory where b is a vector field describing 
the dynamics of the system and / is a running cost. In that framework the 
meaning of v is that of a value function which allows to characterize the 
optimal trajectories (here they play the role of characteristic curves). The 
interested reader can find more details on this interpretation in ED- 

4 The Oren-Nayar model (ON—model) 

The diffuse reflectance ON-model gHi sa sa Eoi is an extension of the 
previous L-model which explicitly allows to handle rough surfaces. The 
idea of this model is to represent a rough surface as an aggregation of V- 
shaped cavities, each with Lambertian reflectance properties (see Fig. [^. 
In [48] and, with more details, in m, Oren and Nayar derive a reflectance 
model for several type of surfaces with different slope-area distributions. In 
this paper we will refer to the model called by the authors the “Qualitative 
Model”, a simpler version obtained by ignoring interreflections (see Section 
4.4 of [IHI for more details). 
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dA 


Figure 3: Facet model for surface patch dA consisting of many V-shaped 
Lambertian cavities. Figure adapted from [34] . 

Assuming that there is a linear relation between the irradiance of the image 
and the image intensity, the Id brightness equation for the ON-model is 
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given by 

Id(x) = cos{Oi){A + B s\Tv{a)i&jv{P)M{(pi,(pr)) (14) 

where A = 1 - 0.5cr^(cr^ + 0.33)“^ (15) 

S = 0.45cr2(cr2 + 0.09)“^ (16) 

= max{0, cos((y^^ — (/p^)}. (17) 


Note that A and B are two non negative constants depending on the statis¬ 
tics of the cavities via the roughness parameter a. We set a G [0,7r/2), 
interpreting a as the slope of the cavities. In this model (see Fig. Q, 9i 
represents the angle between the unit normal to the surface N(x) and the 
light source direction a;, Or stands for the angle between N(x) and the ob¬ 
server direction V, ipi is the angle between the projection of the light source 
direction uj and the xi axis onto the (xi, X 2 )-plane, denotes the angle 
between the projection of the observer direction V and the xi axis onto the 
(xi, X 2 )-plane and the two variables a and /3 are given by 

a — max {0^, Or} and (3 = min {0^, Or} . (18) 

Since the vectors u; and V are fixed and given, their projection on the 
incident plane is obtained considering their first two components over three 
(see Eq. ([2T|)). In this way, the quantity max{0 ,cos((/:p^ — ipi)} is computed 
only once for a whole image. 



Figure 4: Diffuse reflectance for the ON-model. Figure adapted from [34] . 
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We define (see Fig. [^: 


cos{9i) 

_ IS! ^ ■ Vm(x) + W3 

yi + |v^x)|2 

(19) 

COs{9r) 

^ Y + ^3 

v/l + |V?/(x)|2 

(20) 

COs{(fr 

- Ifii) = (cJi,u; 2 ) ■ (vi,V 2 ) = -v 

(21) 

sm{9i) 


(22) 

sm{9r) 


(23) 


where 

5c^(V«(x)) := a /1 + |V«(x)|2 - (-W • V«(x) +^3)2 
5„(V«(x)) := a /1 + |Vtl(x)|2 _ . v«(x) + ^3)2. 


For smooth surfaces, we have a = 0 and in this case the ON-model 
reduces to the L-model. In the particular case a; = V = (0,0,1), or, more 
precisely, when cos{(pr — ^i) ^ 0 (e.g. the case when the unit vectors uj and 

V are perpendicular we get cos((^^ — ipi) = —1) the equation (14) simplifies 
and reduces to a L-model scaled by the coefficient A. This means that the 
model is more general and flexible than the L-model. This happens when 
only one of the two unit vectors is zero or, more in general, when the dot 
product between the normalized projections onto the (xi, X 2 )-plane of cv and 

V is equal to zero. 

Also for this diffuse model we neglect the ambient component, setting 
k]j — 1. As a consequence, in the general equation @ the total light 
intensity /(x) is equal to the diffuse component /^(x) (described by the 
equation (jM])). This is why we write /(x) instead of /d(x) in what follows. 

To deal with this equation one has to compute the min and max opera¬ 
tors which appear in (14) and (18). Hence, we must consider several cases 
described in detail in what follows. For each case we will derive a partial 
differential equation that is always a first order nonlinear HJ equation: 


Case 1: 9i > 9r and {cpr — ^i) ^ [0, f) U (|7r, 2ti\ 

The brightness equation © becomes 

/(x) = cos(0j) (cos{(pr - j (24) 
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and by using the formulas ©-(HI) we arrive to the following HJ equation 
I(x)(yi + |Vrt(x)| 2 ) ^ . Vrt(x) - C03) 

_ • v)gu,(Vtt(x))g^(Vn(x))(-(I> • Vrt(x) +CJ3) ^ ^ 

1/1 + |Vrt(x)p (-V • Vn(x) + ^3) 

where u; := ( 001 ^ 002 ) and v := ('r’i,'r’ 2 )- 


Case 2: 0i < 9r and {(fr — ^i) ^ [0, f) U (^tt, 27 r] 

In this case the brightness equation (0 becomes 

/(x) = cos(6»i) (A+Bsm{9r)^^^^^^cos{ipr -(fi) 

V cosie'j) 


and by using again the formulas (19)-(23) we get 


/(x) (1 + 1 Vm(x)|2) 

+A(uj ■ Vn(x) - L 03 ) a/ 1 + |V«(x)p 
-B{u} ■ v) g^{Vu{x)) gy{Vu{x.)) = 0. 


(26) 


(27) 


Case 3: \/9i,6r and {ipr — ffi) G [|, |7r] 

In this case we have the implication max{0, cos{(fr 
ness equation (14) simplifies in 

/(x) = A cos{9i) 


(pi)} = 0. The bright- 
(28) 


and the HJ equation associated to it becomes consequentially 

/(x)(y/l -h |Vi^(x)p) + A{iv • Vi^(x) - Us) = 0, (29) 

that is equal to the L-model scaled by the coefficient A. 


Case 4: 9i = 9r and (pr = p^i 

This is a particular case when the position of the light source a; coincides 
with the observer direction V but there are not on the vertical axis. This 
choice implies max{0, cos{p)i — p^r)} — 1? then defining 9 9i — 9^ — Oi — (3^ 
the equation (14) simplifies to 

/(x) = cos(0) (A-hSsin(0)^ cos(0)“^) (30) 


and we arrive to the following HJ equation 

(/(x) -B){pi + |Vm(x)|2) + A{u; ■ Vn(x) - CJ3) 

^ (-u; ■ Vu(x)+cv3)‘^ ^ ^ 

yi + |V^x)|2 


(31) 
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Note that this four cases are exactly the same cases reported and analyzed in 
[35] . This is not surprising since the reflectance model used there is always 
the same one proposed by Oren and Nayar. However, here we get differ¬ 
ent HJ equations since we consider an orthographic camera projection and 
cartesian coordinates whereas in [35] the HJ equations are derived in spheri¬ 
cal coordinates under a generalized perspective camera projection. Another 
major difference is that in that paper the light source is close to the camera 
but is not located at the optical center of the camera. 

The vertical light case. 

If u; = (0,0,1), independently of the position of V, the analysis is more 
simple. In fact, the first three cases considered above reduce to a single case 
corresponding to the following simplified PDE for the brightness equation 


(14) 


I(x) = 


A 


(32) 


yi + |v«(x)|2' 

In this way we can put it in the following eikonal type equation, analogous 
to the Lambertian eikonal equation 0 : 


|Vit(x)| = /(x) for X e fi, 


(33) 


where 

Following we write the surface as *S(x, 2 ;) = z — u(x) = 0, for x G fl, 
z G M, and VS(x, z) = (— Vi^(x), 1), so becomes 


(/(x)-S) |V5(x,z)| 

+^(—^^^(x, z) ■ oj) 


+B 


( |V5(xj| •^) |V5(x,z)|=0. 


Defining 

and 


d(x,2;) := VS{x^z)/\\/S{x^z)\ 


c(x, z) := /(x) - B H(d(x, z) • 
using the equivalence 


|V*S(x, 2 ;)|= max {a • V*S(x, 2 ;)} 

aedBs 


(35) 


(36) 

(37) 

(38) 


we get 

max {c(x, z) a • V*S(x, z) — Au; • V*S(x, 2 ;)} = 0. (39) 

aedBs 
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(40) 


Let us define the vector field for the ON-model 

^ _ (c(x,z)ai - Aui,c(x,z)a2 - Auj2) 

and 

f^^{x,z,a,v{x)) := -/^u(x)) + 1. (41) 


Then, introducing the exponential Kruzkov change of variable fiv{'x.) = 1 — 
g-/in(x) already done for the L-model, we can finally write the fixed point 
problem in the new variable v obtaining the 


Oren-Nayar Model 



f /iv(x) = r^^(x. 

v(x),Vv), X G D, 

(42) 

{ v(x) = 0, 

X G dfl. 

where 



:= min {hO^(x,a) 

■ Vv(x) + f^^(x, z, a, v(x))} 


aedBs 



Note that the simple homogeneous Dirichlet boundary condition is due to 
the flat background behind the object but a condition like i^(x) = ^(x) can 
also be considered if necessary. Moreover, the structure is similar to the 
previous Lambertian model although the definition of the vector field and 
of the cost are different. 

In the particular case when cos{(pr — ^i) — 0, the equation (14) simply 
reduces to 

/(x) = Acos{9) (43) 

and, as a consequence, the Dirichlet problem in the variable v is equal to 
(42) with c(x, z) = /(x). 


5 The Phong model (PH—model) 

The PH-model is an empirical model that was developed by Phong m in 
1975. This model introduces a specular component to the brightness func¬ 
tion /(x), representing the diffuse component /d(x) in (Q as the Lambertian 
reflectance model. 

A simple specular model is obtained putting the incidence angle equal to 
the reflection one and a;, N(x) and R(x) belong to the same plane. 

This model describes the specular light component /s'(x) as a power of the 
cosine of the angle between the unit vectors V and R(x) (it is the vector rep¬ 
resenting the reflection of the light a; on the surface). Hence, the brightness 
equation for the PH-model is 

/(x) = A)^(N(x) • u;) + ksmx) • V)", (44) 


15 







where the parameter ol G [1,10] is used to express the specular reflection 
properties of a material and and ks indicate the percentages of diffuse 
and specular components, respectively. Note that the contribution of the 
specular part decreases as the value of a increases. 

Starting to see in details the PH-model in the case of oblique light source 
uj and oblique observer V. 

Assuming that N(x) is the bisector of the angle between oj and R(x) (see 
Fig. [^, we obtain 

N(x)= “ + 


which implies 


|a; + R(x)|| 
R(x) = ||a; + R(x)||N(x) 


— u?. 


(45) 


(46) 


From the parallelogram law, taking into account that a;, R(x) and N(x) are 
unit vectors, we can write ||a; + R(x)|| = 2(N(x) • u?), then we can derive 
the unit vector R(x) as follow: 


R(x) = 2(N(x) • a?)N(x) - u; 
-u; • Vi^(x) + 6^3 


= 2 


V v i + iv«(x)r^ ^n(x)-(...,....3) 

-2w-V«( x) + 2w3^ v7 ^ ^ ^ ^ 

1 + |Vu(x)P j (-V«(x), 1) - ( 0 . 1 , 0 - 2 ,0.3). 


(47) 


N 



Figure 5: Geometry of the Phong reflection model. 


With this definition of the unit vector R(x) we have 


R(x) • V = 

/—2a; • V?i(x) + 2 a; 3 \ 

V 1 + |V«(X)|2 ) 


(—Vit(x) • V + Vs) — w • V. 


(48) 
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Then, setting a = 1, that represents the worst case, equation (44) becomes 


/(x)(l+|Vn(x)|2) 

-koi-'^uix) ■ u + W3)(a/ 1 + |V«(x)|2) (49) 

-2ksi-'Vu{x) ■ u; + uis) (-v • Vit(x) + vs) 

+fc5(wV)(l + |V«(x)|2) = 0, 

to which we add a Dirichlet boundary condition equal to zero, assuming 
that the surface is standind on a flat blackground. Note that the cosine in 
the specular term is usually replaced by zero if R(x) • V < 0 (and in that 
case we get back to the L-model). 

As we have done for the previous models, we write the surface as *S(x, z) = 
2 ) — i/(x) = 0, for X E 14, 2 ) G M, and VS'(x, z) = (—Vi^(x), 1), so ( [4^ will be 
written as 


(I(x) + fcs(wV))lV5(x,z)p 

-/ci 5 (VS'(x, 2 :) • u;)(|VS'(x,z)|) (50) 

-2ks{VS{^, z) ■ w)(V5(x, z)-\) = 0. 

Dividing by |VS'(x, z)\, defining d(x, z) as in ( |3^ and c(x) := I(x) + ks{iv ■ 
V), we get 

c(x)|VS'(x,z)! — fc£)(VS'(x,z) • w) (51) 

—2ks(VS{-K, z) ■ w)(d(x, z) • V) = 0. 

By the equivalence |VS'(x, z)| = max {a • VS'(x, z)} we obtain 

aedBs 

max {c(x) a • z) — fc£)(VS'(x, z) • m) (52) 

aedBs 

—2ks{'VS{x, z) ■ u;)(d(x, z) ■ V)} = 0. 

Let us define the vector field 

b^^(x,a) := qph(^^ ^^ M^^(x,z) (53) 

where 

Q^^{x, z) := 2/c5'W3(d(x, z) • V) + (54) 

and 

M[^{x,z) := (c(x)ai - kocvi 

-2ksLOi{d{x,z)-Y)) ,i = 1,2. (55) 

Let us also define 

f^^{x,z,a,v{x)) :=-^^^^(l-/ru(x)) + l. (56) 
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Again, using the exponential Kruzkov change of variable /i'L’(x) = 1 — 

as done for the previous models, we can finally write the nonlinear fixed point 

problem 


Phong Model 

( iiv{x) = r^^(x, 
\ u(x) = 0, 

'u(x), V'L’), X G ri, 

X G dQ. 

(57) 

where 





:= min {b^^(x,a) 

aedBs 

■ Vu(x) + z, a, u(x))} 



A. Oblique light source and vertical position of the observer. 

In the case of oblique light source oo and vertical observer V = (0,0,1), the 
dot product R(x) • V becomes 


R(x) • V 


—2u; • Vi^(x) + 
l + |Vi^(x)p 

—2u; • Vi^(x) + 6 ^ 3(1 — |Vi^(x)p) 
1 + |Vi^(x)|2 


(58) 


The fixed point problem in 
c(x) 

b^^(x, a) 


V will be equal to ( [57| ) with the following choices 
:= /(x)+a;3/c5, 

:= 2ks{d{i^,z) • Lo) (59) 

_ {c{-x.)ai—kD^i^c{-x)a2—kD^2) 


B. Vertical light source and oblique position of the observer. 

When a; = (0, 0,1) the definition of the vector R(x) reported in (47) becomes 

—2Vi^(x) 2 


R(x) = 


1 + |V?/(x)P’ 1 + |Vn(x)p 


1 


(60) 


and, as a consequence, the dot product R(x) • V with general V is 


R(x) • V = 


-2v • Vrt(x) + ^ 3(1 — |Vri(x)|^) 
1 + |Vri(x)p 


Hence, the fixed point problem in v is equal to (57) with 


c(x) 

b^^(x, a) 


/(x) + v^ks, 

2 /cs(d(x,z) • V) + kD, 
QPn\^,z) (c(x)ai,c(x)a 2 ). 


(61) 
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C. Vertical light source and vertical position of the observer. 

If we choose a; = V = (0,0,1) the equation ( [4^ simplifies in 

I(x)(l + lVt((x)p) - ks{l - |Vt((x)p) (63) 

+ |Vn(x)|2) = 0. 

Working on this equation one can put it in the following eikonal type form, 
which is analogous to the Lambertian eikonal equation 0 : 

|Vi/(x)| = /(x) for X G (64) 

where now 


\ — 

kj) - 

2I+(x)/(x) + k‘]^^Q{x) 

(65) 

) — 


( 

\ 


2f/(x) + fcs) 


/+(x 

) ■■= 

I(x) + ks, 

(66) 

/_(x 

) ■■= 

I(x) - kg, 

(67) 

Q(x 

) ■■= 

+ 8kg + 81{x) kg. 

(68) 


Remark on the eontrol interpretation. The above analysis has shown 
that all the cases corresponding to the models proposed by Oren-Nayar and 
by Phong lead to a stationary Hamilton-Jacobi equation of the same form, 
namely 

/i'L’(x) = min {b(x, a) • V'?;(x) + /(x, 2 ;, a, 'r’(x))} 

aedBs 

where the vector field b and the cost / can vary according to the model 
and to the case. This gives to these models a control theoretical interpre¬ 
tation which can be seen as a generalization of the control interpretation 
for the original Lambertian model (which was related to the minimum time 
problem). In this framework, v is the value function of a rescaled (by the 
Kruzkov change of variable) control problem in which one wants to drive 
the controlled system governed by 

y{t) = h{y{t),a{t)) 

(a(-) here is the control function taking values in dBs) from the initial po¬ 
sition X to the target (the silhouette of the object) minimizing the cost 
associated to the trajectory. The running cost associated to the position 
and the choice of the control will be given by /. More informations on this 
interpretation, which is not crucial to understand the application to the SfS 
problem presented in this paper, can be found in m 
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6 Semi-Lagrangian Approximation 


Now, let US state a general convergence theorem suitable for the class of 
differential operators appearing in the models described in the previous sec¬ 
tions. As we noticed, the unified approach presented in this paper has the 
big advantage to give a unique formulation for the three models in the form 
of a fixed point problem 

= T^(x, 'L’, V'L’), for X E ri, (69) 

where M indicates the model, i.e. M = ON^ PH. 

We will see that the discrete operators of the ON-model and the PH- 
model described in the previous sections satisfy the properties listed here. 
In order to obtain the fully discrete approximation we will adopt the semi- 
Lagrangian approach described in the book by Falcone and Ferretti [2T]. The 
reader can also refer to m for a similar approach to various image processing 
problems (including nonlinear diffusions models and segmentation). 

Let Wi = w{xi) so that W will be the vector solution giving the ap¬ 
proximation of the height u at every node Xi of the grid. Note that in one 
dimension the index i is an integer number, in two dimensions i denotes a 
multi-index, i = (^ 1 ,^ 2 )- We consider a semi-Lagrangian scheme written in 
a fixed point form, so we will write the fully discrete scheme as 

Wi = fi(W). (70) 

Denoting by G the global number of nodes in the grid, the operator cor¬ 
responding to the oblique light source is T : ^ that is defined 

componentwise by 

fi(W)-.= min {e-^^^I[W](xj)-TF(xi,z,a)}+ T ( 71 ) 

aedBs 

where I[W] represents an interpolation operator based on the values at the 
grid nodes and 

xf := Xi + hb{xi,a) ( 72 ) 

T ( 73 ) 

F{xi, z, a) := P{xi, 2 ) 03(1 - nWi) ( 74 ) 

P : D X M ^ M is continuous and nonnegative. (75) 

Since w{xi -\- hb{xi^a)) is approximated via I\W] by interpolation on W 
(which is defined on the grid G), it is important to use a monotone interpo¬ 
lation in order to preserve the properties of the continuous operator T in the 
discretization. To this end, the typical choice is to apply a piecewise linear 
(or bilinear) interpolation operator /i[VF] : D ^ M which allows to define a 
function defined for every x E D (and not only on the nodes) 

w{x)=h[W](x)^Y.^i^(^)^J ( 76 ) 

j 
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where 


(77) 


Xij{a) = 1 for X = Xij{a)xj. 

3 3 

A simple explanation for @-((73 is that the coefficients Xij{a) represent 
the local coordinates of the point x with respect to the grid nodes (see m 
for more details and other choices of interpolation operators). Clearly, in 
0 we will apply the interpolation operator to the point x^ = Xi + hb{xi^ a) 
and we will denote by w the function defined by /i[VC]. 

Comparing 0 with its analogue for the vertical light case we can immedi¬ 
ately note that the former has the additional term z, a) which requires 

analysis. 

Theorem 6.1 Let TiiW) the i-th eomponent of the operator defined as in 
0. Then, the following properties hold true: 

1. Let 

as = arg min {e~^^w{xi + hb{xi, a)) — rF^Xi, z, a)} and assume 

aedBs 


P{xi,z)as<l. (78) 

Then 0 < W < implies 0 < T(W) < 

2. T is a monotone operator, i.e., W <W implies T(W) < T(W). 


3. T is a eontraetion mapping in L^([0, l//i)^) if 
P{xi,z)as < /i. 

Proof 

1. To prove that W < ^ implies T{W) < we just note that 

_ -fih 1 

T{W) < -+T = -. (79) 

ji fi 

Let VL > 0; then 

f{W) > -tP{x^, z) a3(l - f^Wi) + T 
= r (1 - P{xi, z) 03(1 - nWi)) . 

This implies that T(W) > 0 if P{xi,z) 03 < 1 since 0 < 1 — fiWi < 1. 
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2. In order to prove that T is monotone, let us observe first that for each 
couple of functions wi and W 2 such that wi{x) < W 2 {x) for every x G 12 
implies 

e~^^[wi{x + hb{x^ a*)) — W 2 {x + hb{x^a))] (81) 

-rP{x,z){al{l - fiwi{x)) - as{l - /j.W 2 {x))) 

< e~^^[wi{x + hb{x^a)) — W 2 {x + hb{x^a))] 

+tP(x, z)as{wi{x) - W 2 (x)) 

where a* and a are the two arguments corresponding to the minimum 
on a of the expression 

e~^^w{x + hb{x^ a)) — rP{x^ z)as{l — jxw{x)) (82) 

respectively for w = wi^W 2 - Hence, if we take now two vectors W 
and W such that W <W and we denote respectively by w and w the 
corresponding functions defined by interpolation, we will have w{x) — 
Ii\W]{x) < Ii\W]{x) — w{x)^ for every x G 12, because the linear 
interpolation operator Ii is monotone. Then, by setting = 

Xi + hb{xi^a) we get 

fi{W) - %{W) < e-^[h[W]{xt)) - h[W]{xi))] 
+TP{xi,z)az{Wi - Wi) < 0 ( 83 ) 


where the last inequality follows from P > 0. So we can conclude that 
T{W) — T(W) < 0. Note that this property does not require condition 
(1^ to be satisfied. 


3. Let us consider now two vectors W and W (dropping the condition 
W <W) and assuming 


P{xi,z)az < fi. 


(84) 


To prove that T is a contraction mapping note that following the 
same argument used to prove the second statement we can obtain 
(83). Then, by applying the definition of /i, we get 

fiW) - f{W) < (e-^^ + TP{xi, z)a3) 11- W| loo. 

Reversing the role of W and VL, one can also obtain 

f{W) - f{W) < (e-'^^ + rPixi, z)a3) 11- W| U 

and conclude then that T is a contraction mapping in if and only 
if 

^-fih z)a3 < 1 (85) 

and this holds true if the bound (ISdl) is satisfied. 
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□ 

Remark on the role of fi. The parameter fi can be tuned to satisfy the 
inequality which guarantees the contraction map property for T. This pa¬ 
rameter adds a degree of freedom in the Kruzkov change of variable and 
modifies the slope of v. However, in the practical applications we have done 
in our tests, this parameter has been always set to 1 so in our experience 
this parameter does not seem to require a fine tuning. 


Remark on the ehoiee of the interpolation operator. Although Ii can 
be replaced by a high-order interpolation operator (e.g. a cubic local La¬ 
grange interpolation), the monotonicity of the interpolation operator plays 
a crucial role in the proof because we have to guarantee that the extrema 
of interpolation polynomial stay bounded by the minimum and maximum 
of the values at the nodes. This property is not satisfied by quadratic or 
cubic local interpolation operators and the result is that this choice intro¬ 
duces spurious oscillations in the numerical approximation. A cure could 
be to adopt Essentially Non Oscillatory (ENO) interpolations. A detailed 
discussion on this point is contained in m- 

Remark on the minimization. In the definition of the fixed point operator 
T there is a minimization over a G dB^. A simple way to solve it is to 
build a discretization of dB^ based on a finite number of points and get the 
minimum by comparison. One way to do it is to discretize the unit sphere 
by spherical coordinates, even a small number of nodes will be sufficient 
to get convergence. A detailed discussion on other methods to solve the 
minimization problem is contained in m- 

Let us consider now the algorithm based on the fixed-point iteration 

given. 



We can state the following convergence result 


Theorem 6.2 Let be the sequenee generated by (86). 
ing results hold: 

1. LetW^eS = {WeR^:W> f{W)}, then the 


Then the follow- 
eonverges mono- 


tonieally deereasing to a fixed point VE* of the T operator; 

2. Let us ehoose ja > 0 so that the eondition P{xi^z)as < fi is satisfied. 
Then, eonverges to the unique fixed point VE* of the T. Moreover, if 
G S the eonvergenee is monotone deereasing. 


Proof. 
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1. Starting from a point in the set of super-solutions 5, the sequence is 
non increasing and lives in S which is a closed set bounded from below 
(by 0). Then, converges and the limit is necessarily a fixed point 
for f. 


2. The assumptions guaranteed by Theorem 


6.1 


are satisfied and T is 


a contraction mapping in [0,l//i], so the fixed point is unique. The 
monotonicity of T implies that starting from G S the convergence 
is monotone decreasing. 


□ 


It is important to note that the change of variable allows for an easy choice 
of the initial guess G S for which we have the natural choice — 
(l//i, l//i,..., l//i) and monotonicity can be rather useful to accelerate con¬ 
vergence as shown in j20| . A different way to improve convergence is to 
apply Fast Sweeping or Fast Marching methods as illustrated in miEg. A 
crucial role is played by boundary conditions on the boundary of fl, where 
usually we impose the homogeneous Dirichlet boundary condition, 'L’ = 0. 
This condition implies that the shadows must not cross the boundary of fl, 
so the choice = 0 corresponding to an infinite shadow behind the surface 
is not admissible. However, other choices are possible: to impose the height 
of the surface on dO. we can set v — g or to use a more neutral boundary 
condition we can impose v — 1 (state constraint boundary condition). More 
informations on the use of boundary conditions for these type of problems 
can be found in m- 


6.1 Properties of the discrete operators and 


We consider a semi-Lagrangian (SL) discretization of (42) written in a fixed 
point form, so we will write the SL fully discrete seheme for the ON-model 
as 

( 87 ) 


Wi = Tl^^\W), 


where ON is the acronym identifying the ON-model. Using the same nota¬ 
tions of the previous section, the operator corresponding to the oblique light 
source is ton . 

^ that with linear interpolation can be written as 




min 

aedBs 


{e-'^Vi[W](x+) - rF^^ixi, z,a)} + r, 
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where 


T := 


a) := {c{xi, z)ai - Aoji, c{xi, z)a2 - Auj2) 
Aids 

c{xi, z) := I(xi) - B + B(d(xi, z) • u;)^ 
d(xi,z) := VS(xi,z)/lVS(xi,z)\ 

F^^(xi, z, a) := P^^{xi, z)a3{l - fiWi)) 


1 - 

fi 


jON 


(xi,z) := 


c{xi,z) 

Aujs 


Note that, in general, will not be positive but that condition can be 

obtained tuning the parameter a since the coefficients A and B depend on 
a. This explains why in some tests we will not be able to get convergence for 
every value of a G [0,7r/2). Once the non negativity of is guaranteed, 
we can follow the same arguments of Theorem 6T to check that the discrete 
operator satisfies the thre e properties which are necessary to guarantee 
convergence as in Theorem 6 ^ provided we set P = pON statement. 

For the Phong model, the semi-Lagrangian discretization of (57) written 
in a fixed point form gives 

Wi = fP(W), ( 88 ) 

where that is defined componentwise by 

fP{W) min {e-^'^Ii[W]{x+) - {xi, z,a)} + r, 

aedBs 

where, in the case of oblique light source and vertical position of the observer, 

1 




b^^{xi,a) := 




{c{xi)ai - kDC0i,c{xi)a2 - kDC02) 


Q^^{xi,z) 

c{xi) := I{xi) Fcosks 
d{xi,z) := VS{xi,z)/\VS{xi,z)\ 

Q^^{xi, z) := 2 ks{d{xi, z) ■ u) + kouJz 
F^^{xi, z, a) := P^^{xi, 2 ) 03(1 - nWi)) 

Here the model has less parameters and P^^ will always be nonnegative. 
Again, following the same arguments of Theorem |6.1[ we can check that the 
discrete operator satisfies the t hree properties which are necessary to 
guarantee convergence as in Theorem 6.2 provided we set P = P^^ in that 
statement. 
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7 Numerical Simulations 


In this section we show some numerical experiments on synthetic and real 
images in order to analyze the behavior of the parameters involved in the 
ON-model and the PH-model and to compare the performances of these 
models with respect to the classical L-model and with other numerical meth¬ 
ods too. All the numerical tests in this section have been implemented in 
language C++. The computer used for the simulations is a MacBook Pro 
13” Intel Core 2 Duo with speed of 2.66 GHz and 4GB of RAM (so the CPU 
times in the tables refer to this specific architecture). 

We denote by Q the discrete grid in the plane getting back to the double 
index notation G := card{Q) = n x m. We define Gin ~ {xij : Xij G 
D} as the set of grid points inside D; Gout := G \ Gin- The boundary 

will be then approximated by the nodes such that at least one of the 
neighboring points belongs to Gin- For each image we define a map, called 
mask^ representing the pixels Xij G Gin in white and the pixels Xij G Gout 
in black. In this way it is easy to distinguish the nodes that we have to use 
for the reconstruction (the nodes inside D) and the nodes on the boundary 

(see e.g. Fig. |6(b)[ ). 

Regarding the minimization over a G dB^ that appears in the definition 
of the fixed point operators associated to the models, in all the tests we 
discretize the unit sphere by spherical coordinates, considering 12 steps in 
9 and 8 in 0, where 9 is the zenith angle and (j) is the azimuth angle. 

7.1 Synthetic tests 

If not otherwise specified, all the synthetic images are defined on the same 
rectangular domain containing the support of the image, D = [—1,1] x 
[—1,1]. We can easily modify the number of the pixels choosing different 
values for the steps in space Ax and Ay. The size used for the synthetic 
images is 256 x 256 pixels, unless otherwise specified. X and Y represent 
the real size (e.g. for Q = [—1,1] x [—1,1], X = 2, Y = 2). For all the 
synthetic tests, since we know the algebraic expression of the surfaces, the 
input image rendering in gray-levels is obtained using the corresponding 
reflectance model. This means that for each model and each value of pa¬ 
rameter involved in it, the reconstruction will start from a different input 
image. Clearly, this is not possible for real images, so for these tests the 
input image will be always the same for all the models, independently of the 
values of the parameters. Moreover, we fix /i = 1 and we choose the value 
of the tolerance y for the iterative process equal to 10“^ for the tests on 
synthetic images, using as stopping rule — W^\max ^ 9^ where A: + 1 

denotes the current iteration. We will see that dealing with real images, 
sometimes we will need to increase 77. 
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Test 1: Sphere. For this first test we will use the semisphere defined as 


u{x,y) = - y 2 (^x,y)eGin, 

u{x,y) = 0 {x,y)eGout, 

where 

Gin ■■= {{x,y) :x^ + y^ < r^}, 

r:=2i2f^+2S. 

Zi 

and 5 := max{Ax, Ay}. 

As example, we can see in Fig. the input image, the corresponding 
mask and the surface reconstructed by the L-model. The values of the 

on 

(a) Sphere Input (b) Sphere Mask 

Figure 6: Sphere via the L-model: (a) Input image; (b) Mask; (c) surface. 

parameters used in the simulations are indicated in Table Note (in Table 
1^ that when the specular component is zero for the PH-model, we just have 
the contribution of the diffuse component so we have exactly the same error 
values of the L-model, as expected. By increasing the value of the coefficient 
ks and, as a consequence, decreasing the value of ko: in the PH-model the 
LP‘{I) and L^{I) errors on the image grow albeit slightly and still remain of 
the same order of magnitude, whereas the errors on the surface decrease. For 
the ON-model, the same phenomenon appears when we set the roughness 
parameter a to zero: we bring back to the L-model and, hence, we obtain 
the same errors on the image and the surface. Note that the errors in L?{I) 
and L^{I) norm for the image and the surface, decrease by increasing the 
value of a. This seem to imply that the model and the approximation work 
better for increasing roughness values. 

We point out that the errors computed on the images / in the different 
norms are over integer between the input image and the image computed a 
posteriori using the value of u just obtained by the methods. This is why 
small errors becomes bigger since can jump from an integer to the other one, 
as for the Phong case. 

In Table we reported the number of iterations and the CPU time (in 
seconds) referred to the three models with the parameter indicated in Table 



(c) Sphere surface 


(90) 

(91) 


27 





Table 1: Sphere: parameter values used in the models. 


Model 

a 

ko 

ks 

a 

LAM 

ON-00 

0 




ON-04 

0.4 




ON-08 

PH-sOO 

0.8 

1 

0 

1 

PH-S04 


0.6 

0.4 

1 

PH-sOS 


0.2 

0.8 

1 


For all the models, also varying the parameters involved, the number 
of iterations is always about 2000 e the CPU time slightly greater than 2 
seconds, so the computation is really fast. 


Table 2: Synthetic sphere: iterations and CPU time in seconds for the 
models with vertical light source u; = (0,0,1). 


SL-Schemes 

Iter. 

[sec.] 

LAM 

2001 

2.14 

ON-00 

2001 

2.06 

ON-04 

2020 

2.24 

ON-08 

2016 

2.24 

PH-sOO 

2001 

2.11 

PH-S04 

2008 

2.45 

PH-sOS 

2056 

2.27 


Table 3: Synthetic sphere: L^, errors with vertical light source a; = 
( 0 , 0 , 1 ). 


SL-Schemes 

L2(/) 


L\S) 

L°°(5) 

LAM 

0.0046 

0.0431 

0.0529 

0.0910 

ON-00 

0.0046 

0.0431 

0.0529 

0.0910 

ON-04 

0.0039 

0.0353 

0.0513 

0.0882 

ON-08 

0.0035 

0.0314 

0.0506 

0.0881 

PH-sOO 

0.0046 

0.0431 

0.0529 

0.0910 

PH-S04 

0.0064 

0.0471 

0.0511 

0.0896 

PH-sOS 

0.0090 

0.0706 

0.0386 

0.0752 


Test 2: Ridge tent. In tests on synthetic images, the relevance of the 
choice of a model depends on which model was used to compute the images. 
In the previous test, the parameters that are used for the 3D reconstruction 
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are identical to those used to compute the synthetic sphere input images, 
so there is a perfect match. However, for real applications, it is relevant to 
examine the influence of an error in the parameter values. To this end we 
can produce an input image with the Oren-Nayar model using a = 0.1 and 
then process this image with the same model using a different value of a to 
see how the results are affected by this error. This is what we are going to 
do for the ridge tent. Let us consider the ridge tent defined by the following 
equation 

{u{x,y) = min{-2|a;| + -|y| + 

< {x,y)eGin, (92) 

^u{x,y) = 0 ix,y)eGout, 

where 

In Fig. [^we can see an example of reconstruction obtained by using the ON- 
model with a = 0.3, under a vertical light source a; = (0,0,1). A first remark 


o 

(a) Tent Input 

Figure 7: Tent via the ON-model with a = 0.3: (a) Input image; (b) 3D 
reconstruction. 

is that the surface reconstruction is good even if in this case the surface is not 
differentiable. Moreover, note that there are no oscillations near the kinks 
where there is jump in the gradient direction. Let us examine the stability 
with respect to the parameters. We have produced seven input images for 
the ridge tent, all of size 256 x 256, with the following combinations of models 
and parameters: 

LAM Lambertian model; 

ONI Oren-Nayar model with cr = 0.1; 

ON3 Oren-Nayar model with a = 0.3; 

ON5 Oren-Nayar model with a = 0.5; 



(b) Tent surface 
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PHI Phong model with a = 1 and = 0.1; 

PH3 Phong model with a = 1 and ks = 0.3; 

PH5 Phong model with a = 1 and ks = 0.5. 

Then we have computed the surfaces corresponding to all the parameter 

choices (i.e. matching and not matching the first choice). The results ob¬ 
tained in this way have been compared in terms of L? and norm errors 
with respect to the original surface. The errors obtained by the ON-model 
are shown in Tableland in Table [5] for the PH-model. 


Table 4: Ridge tent, ON-model: errors for the surface. In each 

column the model used to produce the input image, in the row the model 
used for the reconstruction. 


L2 

LAM 

ONI 

ON3 

ON5 

LAM 

0.0067 

0.0172 

0.0933 

0.1920 

ONI 

0.0082 

0.0068 

0.0821 

0.1801 

ON3 

0.0832 

0.0700 

0.0086 

0.1033 

ON5 

0.1946 

0.1784 

0.0923 

0.0067 

LOO 

LAM 

ONI 

ON3 

ON5 

LAM 

0.0094 

0.0315 

0.1942 

0.4060 

ONI 

0.0199 

0.0093 

0.1701 

0.3805 

ON3 

0.1784 

0.1507 

0.0118 

0.2156 

ON5 

0.4104 

0.3769 

0.1976 

0.0094 


Table 5: Ridge tent, PH-model: errors for the surface. In each 

column the model used to produce the input image, on the row the model 
used for the reconstruction. 


L2 

LAM 

PHI 

PH3 

PH5 

LAM 

0.0067 

0.0841 

0.2867 

0.6146 

PHI 

0.0586 

0.0067 

0.1996 

0.5031 

PH3 

0.1403 

0.0955 

0.0073 

0.2687 

PH5 

0.1915 

0.1664 

0.0976 

0.0060 

LOO 

LAM 

PHI 

PH3 

PH5 

LAM 

0.0094 

0.1740 

0.6068 

1.3123 

PHI 

0.1243 

0.0108 

0.4202 

1.0741 

PH3 

0.3167 

0.2245 

0.0149 

0.5718 

PH5 

0.4503 

0.4241 

0.2907 

0.0093 
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Analyzing the errors in Table and Table we can observe that using 
the same model to generate the input image and to reconstruct the surface 
is clearly the optimal choice. The errors on the surface grow more as we 
consider a parameter a other than the one used to generate the input image 
as data for the 3D reconstruction. For the ON-model we loose one or two 
order of magnitude, depending on the “distance” of the parameter from the 
source model. For the PH-model we can observe that the L? and errors 
grow more as we consider a different ks for the generation of the image and 
for the reconstruction, loosing one or two order of magnitude. However, 
the two models seems to be rather stable with respect to a variation of 
the parameters since the errors do not increase dramatically varying the 
parameters. 

Test 3: Concave/convex ambiguity for the ON—model. We consider 
this test in order to show that the ON-model is not able to overcome the 
concave/convex ambiguity typical of the SfS problem although it is a model 
more realistic than the classical L-model. Let consider the following function 

u{x,y) = l if(x2 + y2)<2, (93) 

I 0 otherwise. 



Figure 8 : Example of concave/convex ambiguity for the ON-model with 
(j — 0.5 and u; = (0,0,1): (a) Original Surface; (b) Maximal Solution; 
(c) Approximated surface with value in the origin equal to zero. 


We discretize the domain D = [—1.5,1.5] x [—1.5,1.5] with 151 x 151 nodes. 
The fixed point has been computed with an accuracy oi rj — 10“^ and the 
stopping rule for the algorithm based on the fixed point iteration defined in 
mis\W^+^-W^\max < 77 , where k-\-l denotes the current iteration. The 
iterative process starts with = 0 on the boundary and — 1 inside in 
order to proceed from the boundary to the internal constructing a monotone 
sequence (see [TIES] for details on the approximation of maximal solutions). 
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Looking at Fig. [^we can note that the scheme chooses the maximal viscosity 
solution stopping after 105 iterations, which does not coincide with the 
original surface. In order to obtain a reconstruction closer to the original 
surface, we fix the value in the origin at zero. In this way we forced the 
scheme to converge to a solution different from the maximal solution (see 
Fig. |8(c)| ). We obtained this different solution shown in Fig. |8(c)| after 82 
iterations. 

Test 4: Concave/convex ambiguity for the PH—model. The fourth 
synthetic numerical experiment is related to the sinusoidal function defined 
as follow 

u{x, y) = 0.5 + 0.5 sm(^) sm(^), 

^ {x,y)eGin, 

^u{x,y) = 0, ix,y)eGout- 

With this test we want to show that also the PH-model is not able to 
overcome the concave/convex ambiguity typical of the SfS problem. 


in out reconstruction 




Figure 9: Synthetic sinusoidal function: example of concave/convex ambi¬ 
guity for the PH-model with ks = 0,0.5, 0.8 for each row from the top to 
the bottom, respectively. 


Fig. H shows the results related to the PH-model with ks = 0,0.5, 0.8. 
In the first column one can see the input images generated by the PH-model 
using the values of the parameter before mentioned. In the second column 
we can see the output images computed a posteriori using the depth just 
computed and approximating the gradient via finite difference solver. What 
we can note is that even if the reconstructed a posteriori images match 
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with the corresponding input images, the SL method always chooses the 
maximal solution even varying the parameters and ks- By adding some 
informations as shown in the previous test it is possible to achieve a better 
result, but these additional informations are not available for real images. 

Test 5: Role of the boundary conditions. With this fifth test we want 
to point out the role of the boundary condition (BC), showing how good BC 
can significantly improve the results on the 3D reconstruction. We will use 
the synthetic vase defined as follow 

u{x,y) = {x,y) e Gin 

uix, y) = g{x, y) {x, y) G Gout, 

where y := y/Y^ 

P{y) := (-10.8 f + 7.2 f + 6.6 f - 3.8 f 

-1.375 + 0.5^ +0.25) X 


(95) 

(96) 


and 


Gin\^ {{x,y)\P{y) >x}. 


In Fig. 10 one can see on the first row the input images (size 256 x 256) 
generated by the L-model, ON-model and the PH-model, from left to right 
respectively. On the second row we reported the 3D reconstruction with 
homogeneous Dirichlet BC {g{x^y) =0). As we can see, there is a con¬ 
cave/convex ambiguity in the reconstruction of the surface. If we consider 
the correct boundary condition, that is the height of the surface at the 
boundary of the silhouette that we can easily derive in this case being the 
object a solid of rotation, what we obtain is visible in the third row of the 
same Fig. 

In Tables and 0 we can see the number of iterations, the CPU time 
and the error measures in L? and norm for the method used with homo¬ 
geneous and not homogeneous Dirichlet BC respectively. Looking at these 
errors we can note that in each table the values are almost the same for the 
different models. Comparing the values of the two tables one can see that 
we earn an order of magnitude using good BC, that confirms what we noted 
looking at Fig. 


Test 6: Comparison with other numerical approximations. In this 
sixth and last test of this section dedicated to synthetic tests we will compare 
the performance of our semi-Lagrangian approach with other methods used 
in the literature. For this reason, we will use a very common image used in 


the literature, that is the vase, defined in the previous test through the (95). 


More in details, we will compare the performance of our semi-Lagrangian 
method with the Lax-Friedrichs Sweeping (LFS) scheme adopted by Ahmed 
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L—model 


ON—model PH—model 



Figure 10: Synthetic vase under vertical light source (a; = (0,0,1)): ex¬ 
ample of concave/convex ambiguity solved by using correct Dirichlet BC. 
On the first row from left to right: input images generated by L-model, 
ON-model with a — 0.2 and PH-model with ks — 0.4. On the second row: 
3D reconstruction with homogeneous Dirichlet BC. On the last row: 3D 
reconstruction with Dirichlet BC u{x^y) — g{x^y). 


Table 6: Synthetic Vase: Number of iterations, CPU time in seconds and 
L\ LOO qytots on the surface with vertical light source u; = (0,0,1) using 
homogeneous Drichlet BC. 


SL-Schemes 

Iter. 

[sec.] 

lHs) 


LAM 

1010 

0.54 

0.1614 

0.3015 

ON-02 

1011 

0.53 

0.1613 

0.3015 

ON-04 

1008 

0.54 

0.1610 

0.3015 

PH-S02 

1007 

0.54 

0.1619 

0.3015 

PH-S04 

1009 

0.55 

0.1621 

0.3015 
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Table 7: Synthetic Vase: Number of iterations, CPU time in seconds and 
errors on the surface with vertical light source a; = (0,0,1) using 
not homogeneous Drichlet BC. 


SL-Schemes 

Iter. 

[sec.] 

LHS) 

L°°(5) 

LAM 

1337 

0.70 

0.0286 

0.0569 

ON-02 

1335 

0.70 

0.0284 

0.0558 

ON-04 

1341 

0.70 

0.0282 

0.0562 

PH-S02 

1330 

0.70 

0.0284 

0.0560 

PH-S04 

1330 

0.70 

0.0280 

0.0558 


and Farag [T] under vertical and oblique light source. Also these authors 
derive some similar HJ equations for the L-model in [2], and generalize this 
approach for various image conditions in [1], comparing their results on the 
only L-model with the results shown in [63] and the algorithms reported in 
[85] . Unfortunately, we cannot compare our semi-Lagrangian approximation 
for the PH-model with no other schemes since, to our knowledge, there are 
no table of errors for the PH-model under orthographic projection in the 
literature. In order to do the comparison, we will consider the vase image 
of size 128 x 128 as used in the other papers. We start to remind the error 
estimations used: given a vector A representing the reference depth map on 
the grid and a vector A representing its approximation, we define the error 
vector as e = A — A and 


erri := ||e||^i 


err2 := \ \e\\L2 


= -E 

N ^ 




1/2 


where N is the total number of grid points used for the computation, i.e. the 
grid points belonging to Gin- These estimators are called mean and standard 
deviation of the absolute error. In Table we compare the error measures 
for the different SfS algorithms under the L-model with a; = (0,0,1). What 
we can note is that our semi-Lagrangian method is better than the other 
ones also if we consider Dirichlet boundary condition equal to zero (as used 
by Ahmed and Farag in their work [1]), but the better result is with correct 
BC, shown in bold in the last row. In Table we can see the same methods 
applied to the L-model but under a different light source, that is (1,0,1). 
Also in this case, our approach obtains always the smallest errors and the 
best is with not homogeneous Dirichlet BC, as noted before. 

Only with respect to the LFS used by Ahmed and Farag, we can com¬ 
pare the performance of our semi-Lagrangian scheme under the ON-model 
as well, since the other authors only consider the L-model. In this con¬ 


text, we show in Table 10 the error measures for the two SfS algorithms 
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Table 8: Synthetic vase: error measures related to the different methods 
for the L-model with vertical light source oo = (0,0,1). In bold the best 
performances. 


Methods 

erri 

err2 

Best US] 

8.1 

11.1 

m 

2.8 

2.0 

m 

0.22 

0.4 

Our proposed {BC = 0) 

0.1570 

0.1717 

Our proposed {BC 7 ^ 0 ) 

0.0349 

0.0385 


Table 9: Synthetic vase: error measures related to the different methods 
for the L-model with oblique light source (1,0,1). In bold the best perfor¬ 
mances, 


Methods 

erri 

err2 

Best US] 

7.9 

13.9 

m 

4.1 

2.6 

m 

1.2 

2.2 

Our proposed {BC = 0) 

0.0683 

0.1061 

Our proposed {BC 7 ^ 0 ) 

0.0218 

0.0242 

with a = 0.2, under vertical position of light source and viewer (u; = 

(0,0,1),V = (0,0,1)). As 

before, the 

best result is obtained using the 

Table 10: Synthetic vase: error measures 

related to the SL and LFS methods 

for the ON-model with a = 

0.2 under vertical light source (u; = (0,0,1)) 

and vertical viewer V = (0,0,1). In bold the best performances. 

Methods 

erri 

err2 

lU 

0.6 

1.0 

Our proposed {BC = 0) 

0.1568 

0.1715 

Our proposed {BC 7 ^ 0 ) 

0.0348 

0.0384 


semi-Lagrangian scheme with not homogeneous Dirichlet BC. The recon¬ 
structions corresponding to the error measures shown in the three last Ta¬ 
bles [^[9, 10 obtained applying our scheme and compared to [1] are shown in 
Fig. The reconstruction obtained by the two methods are comparable. 


In particular, in the second column regarding the oblique light source case 
we can note that our scheme reconstructs a surface that incorporates the 
black shadow part (see [23] for more details on this technique), avoiding 
the effects of “dent” present in the reconstruction obtained by Ahmed and 
Farag visible in the last row, second column. 

Finally, in Table HH we reported the number of iterations and the CPU 
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L—model L—model ON—model 

vertical oblique vertical 



Figure 11: Synthetic vase: From top to bottom, Input images, recovered 
shapes by our approach with homogeneous Dirichlet BC and with not homo¬ 
geneous BC, recovered shape by [I]. First column: L-model with vertical 
light source (0,0,1). Second column: L-model with oblique light source 
(1,0,1). Third column: ON-model with a — 0.2, a; = (0,0,1), V = (0,0,1). 
Input images size: 128 x 128. 


time in seconds with the comparison with respect to [T]. This shows that the 
SL-scheme is competitive also in terms of CPU time. Of course, in the case 
of oblique light source the number of iterations, and hence the CPU time 
needed is much more bigger. Just think that for the reconstruction under 
the L-model with oblique light source (1,0,1) and BC ^ 0 visible in the 
second column of the third row of Fig. 11, we need 15513 iterations (CPU 
time: 147.40 seconds). If we double the size of the input image, considering 
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Table 11: Synthetic vase: iterations and CPU time in seconds for the L- 
model and the ON-model with vertical light source a; = (0,0,1). Image 
size: 128 x 128. 


Schemes 

Model 

Iter. 

[sec.] 

in 

L-model 

- 

0.5 

Our approach {BC = 0) 

L-model 

611 

0.09 

Our approach {BC ^ 0) 

L-model 

792 

0.11 

m 

ON-model 

- 

1.5 

Our approach {BC = 0) 

ON-model 

612 

0.09 

Our approach {BC ^ 0) 

ON-model 

791 

0.12 


the vase 256 x 256 visible in Fig. for the reconstruction visible in the 
third row, first column of the same Fig. W, we need 30458 iterations to get 
convergence, that we obtain in 1163 seconds. 

Since it is difficult to compare the performance of our approach based 
on the PH-model with other schemes based on the same reflectance model 
under orthographic projection for the consideration made before, we report 
in Fig. [^the performances obtained by our method based on the PH-model 
compared to the approximate Ward’s method on the vase test (Cf. Fig. 7 in 
[1]). What we can see is that both the two methods reconstruct the surface 
in a quite good way, without particular distinctions in the goodness of the 
reconstructions. In order to analyze the results not only in a qualitative 
way but also in a quantitative one, we report the mean and the standard 
deviation of the absolute errors in Table Looking at this table, we can 
note that our approach seems to be superior, obtaining the reconstruction 
with errors smaller of one order with respect to the other model. 


Table 12: Synthetic vase: error measures related to the cases shown in Fig. 


12 In bold the best performances. 


Methods 

erri 

err2 

Ward in [I] 

0.8 

1.3 

PH—model 

0.03 

0.04 


7.2 Real tests 


In this subsection we consider real input images. We start considering the 
golden mask of Agamemnon taken from 17H1 and then modified in order to 
get a picture in gray tones. The size of the modified image really used is 
507 X 512. The input image is visible in Fig. 13(a)[ the associated mask used 
for the 3D reconstruction in Fig. 13(b)[ The second real test is concerning 
the real vase (RV in the following) visible in Fig. 16, taken from [19]. The 
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reconstruction 


Input 



Figure 12: Synthetic vase: The first row shows the input image and the 
recovered shape by our approach with not homogeneous BC based on the 
PH-model with ks = 0.2, a; = (0,0,1) and V = (0,0,1). The second row 
shows the input image and the recovered shape by the approximate Ward’s 
method illustrated in [1] with a = 0.2, pd — 0.67, ps = 0.075, oo — (0,0,1), 
V = (0,0,1). Input images size: 128 x 128. 


size of the input image shown in Fig. |16(a)| is 256 x 256. The reconstruction 
domain shown in Fig. |16(b)| is constituted by the pixels situated on 

the vase. For the real cases, the input image is the same for all the models 
and we can compute only errors on the images since we do not know the 
height of the original surface. For the real tests we will use the same stopping 
criterion for the iterative method before defined for the synthetic tests, i.e. 

< V- 


Test 7: Agamennon mask. For this test we will compare the results 
regarding 3D reconstruction of the surface obtained with a vertical light 
source ojyert = (0,0,1) and an oblique light source cjobl — (0,0.0995,0.9950). 
The values of the parameters used in this test are reported in Table For 


a vertical light source, we refer to Table 14 for the number of iterations and 


the CPU time (in seconds) and to Table 15 for the errors obtained with a 
tolerance p = 10“^ for the stopping rule of the iterative process. Clearly, the 
number of iteration and the errors of the two non-Lambertian models are 
the same of the classical Lambertian model when a for the ON-model and 
ks for the PH-model are equal to zero (and for this reason we do not report 
them in the tables). In all the other cases, the non-Lambertian models are 
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(a) Agamemnon Input (b) Agamemnon Mask 

Figure 13: Agamemnon images (size 507 x 512): (a) Input image; (b) Mask. 

Table 13: Real Agamemnon mask: parameter values used in the models 
with vertical light source a; = (0,0,1). 


Model 

a 

ko 

ks 

a 

LAM 

ON-04 

0.4 




ON-08 

0.8 




ON-10 

PH-S04 

1 

0.6 

0.4 

1 

PH-sOS 


0.2 

0.8 

1 

PH-slO 


0 

1 

1 


faster in terms of CPU time and need a lower number of iterations with 
respect to the L-model. 


Table 14: Real Agamemnon mask: iterations and CPU time in seconds for 
the models with vertical light source oo = (0, 0,1). 


SL-Schemes 

Iter. 

[sec.] 

LAM 

3921 

24.48 

ON-04 

2751 

12.48 

ON-08 

1943 

11.41 

ON-10 

1818 

8.79 

PH-S04 

2127 

9.89 

PH-sOS 

1476 

6.90 

PH-slO 

1325 

6.33 


In Table ^ we can observe that the errors produced by the ON-model 
increase by increasing the value of a. However, the errors are lower than 
the error obtained with the Lambertian model. With respect to the PH- 
model, all the errors increase by increasing the value of the parameter ks^ 
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as observed for synthetic images. 


Table 15: Real Agamemnon mask: errors with vertical light source 


SL-Schemes 

L\I) 

L^{I) 

LAM 

0.0371 

0.4745 

ON-04 

0.0375 

0.4627 

ON-08 

0.0440 

0.4627 

ON-10 

0.0501 

0.4627 

PH-S04 

0.0383 

0.4824 

PH-sOS 

0.0391 

0.4941 

PH-slO 

0.0393 

0.5098 


In Fig. [T^we can see the output image and the 3D reconstruction in a 
single case for each models. What we can note is that no big improvements 
we can obtain visually. 


LAM ON4 PH4 



Figure 14: Agamennon mask: results with vertical light source. On the 
first row the output images, on the second row the 3D reconstruction with 
vertical view, on the third row the 3D reconstruction with oblique view. 
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For the oblique light case, we consider the values for the parameters 


reported in Table 16 


Table 16: Real Agamemnon mask: parameter values used in the models 
with an oblique light source (jJoM = (0,0.0995,0.9950). 


Model 

a 

ko 

ks 

a 

LAM 

ON-01 

0.1 




ON-02 

0.2 




ON-03 

PH-S02 

0.3 

0.8 

0.2 

1 

PH-S03 


0.7 

0.3 

1 

PH-S04 


0.6 

0.4 

1 


Looking at Table we can note that the oblique cases require higher 
CPU time with respect to the vertical cases due to the fact that the equations 
are more complex because of additional terms involved. Because of these 
additional terms involved in the oblique case, in Table we have reported 
the results obtained using the parameters shown in Table W with a value of 
the tolerance r] for the stopping rule of the iterative method equal to 10“^. 
This is the maximum accuracy achieved by the non-Lambertian models since 
roundoff errors coming from several terms occur and limit the accuracy. 


Table 17: Real Agamemnon mask: number of iterations and CPU time 
in seconds for the different models with oblique light source ojobl = 
(0,0.0995,0.9950). 


SL-Schemes 

Iter. 

[sec.] 

LAM 

321 

117.9 

ON-01 

315 

246.0 

ON-02 

361 

281.5 

ON-03 

396 

264.6 

PH-S02 

427 

285.2 

PH-S03 

564 

373.6 

PH-S04 

680 

484.1 


In Fig. we can see the output image and the 3D reconstruction in 
a single case for each models. What we can note is that also using more 
realistic illumination models as the two non-Lambertian considered, we do 
not obtain a so better approximation of the ground truth solution. This is 
due to the missing important informations (e.g. the correct oblique light 
source direction, the values of the parameter involved that are known for 
the models but not available in the real cases) and also due to the fact that 
we are using a first order method of approximation. 
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Table 18: Real Agamemnon mask: L?^ errors with oblique light source 
= (0,0.0995,0.9950). 


SL-Schemes 

L\I) 

L°°{I) 

LAM 

0.0585 

0.4863 

ON-01 

0.0663 

0.4588 

ON-02 

0.0670 

0.4471 

ON-03 

0.0708 

0.5451 

PH-S02 

0.1141 

0.5725 

PH-sOS 

0.1580 

0.6196 

PH-S04 

0.2063 

0.6706 


LAM ON4 PH4 



Figure 15: Agamennon mask: results with oblique light source ujobl — 
(0,0.0995,0.9950). On the first row the output images, on the second row 
the 3D reconstruction with vertical view, on the third row the 3D recon¬ 
struction with oblique view. 


Test 8: Real Vase. With this test we want to investigate the stability 
of our method with respect to the presence of noise. In fact, looking at the 
Fig. |16(a)| we can consider that RV is a noisy version of the synthetic vase 
used in the Tests 5 and 6. The test was performed using a vertical light 
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(a) Real Vase Input 


(b) Real Vase Mask 


Figure 16: Real Vase images (size 256 x 256): (a) Input image; (b) Mask. 


source uj = (0,0,1). The output images, computed a posteriori by using 
the gradient of u approximated via centered finite differences starting from 
the values of u just computed by the numerical scheme, are visible in Fig. 
HH first column. The reconstruction obtained with the three models are 
visible in the same Fig. 


17, second column. What we can note is that all 


the reconstructions, obtained using homogeneous Dirichlet BC, suffer for a 
concave/convex ambiguity, as already noted for the synthetic vase (SV in 
the following). Since it is visible looking at the SV and the RV tests that 
we obtain results with errors of the same order of magnitude around 10“^, 
considering that RV is a noisy version of SV, this shows that the method is 
stable in the presence of noise in the image. 


Table 19: Real Vase: L^, errors with vertical light source u; = (0, 0,1). 


SL-Schemes 

L\I) 


LAM 

0.0093 

0.0784 

ON-02 

0.0094 

0.0784 

ON-04 

0.0111 

0.0824 

PH-S02 

0.0095 

0.0824 

PH-S04 

0.0098 

0.0824 


Finally, in Fig. 18 one can note the behavior of the ON-model by varying 
the value of the parameter a. Since in real situations we do not know it (as 
other parameters like the light source direction) we can only vary it in order 
to see the one which gives the best fit with the image. Looking at Fig. 
what we can observe is that increasing the value of a the reconstruction 
shows a wider concave/convex ambiguity, which affects more pixels. But 
this holds only in this specific case, not in general for all the real images as 
a rule. 
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Out 3D reconstruction 



Figure 17: Real Vase: Output images and 3D reconstructions. On the first 
row the L-model, on the second row the ON-model with a — 0.2, on the 
third row the PH-model with ks — 0.2. 


ON cr = 0.2 ON cr = 0.4 ON cr = 0.6 



Figure 18: Real Vase: 3D reconstructions related to the ON-model, varying 
a. From left to right a — 0.2,0.4, 0.6. 


Test 9: Corridor. As illustrative example, let us consider a real image of 
a corridor (see Fig. 19) as a typical example of a scene which can be useful for 
a robot navigation problem. The test has been added as illustrative example 
to show that even for a real scene which does not satisfies all the assumptions 
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and for which several informations are missing (e.g. boundary conditions) 
the method is able to compute a reasonable accurate reconstruction in the 
central part of the corridor (clearly the boundaries are wrong due to a lack 
of information). The size of this image is 600 x 383. Note that for this 
picture we do not know the parameters and the light direction in the scene. 
It seems that there is a diffused light and more than one light source. So this 
picture does not satisfy many assumptions we used in the theoretical part. 
In order to apply our numerical scheme we considered a Dirichlet boundary 
condition equal to zero at the wall located at the bottom of the corridor. In 
this way we have a better perception of the depth of the scene. In Fig. [^we 


can see the output images (on the first column) and the 3D reconstructions 
(on the second column) obtained by L-model, ON-model with a = 0.1 and 
PH-model with = 0.2. In this example the PH-model seems to recognize 
the scene better than the ON-model. In some sense this is probably due to 
the fact that it has less parameters so it is easier to tune to a real situation 
where the information on the parameters is not available. We point out that 
this test is just an illustration of the fact that coupling SfS with additional 
informations (e.g. coming from distance sensors to fix boundary conditions) 
can be useful to describe a scene. 



Figure 19: Image of a real scene (size 600 x 383). 


Test 10: Other tests on real images. In order to demonstrate the 
applicability of our proposed approach for real data, we added here other 
experiments conducted on a real urn, real rabbit and real Beethoven’s bust. 
The input images and the related recovered shapes obtained by the three 
models studied under different light directions and parameters are shown 
in Fig. As visible from Fig. the results are quite good, even for 


pictures like the rabbit or the bust of Beethoven, which have many details, 
even if they still suffer for the concave/convex ambiguity typical of the SfS 
problem. 
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Out 


3D reconstruction 



Figure 20: Output images and 3D reconstructions of a scene for a robot 
path planning application. On the first row the L-model, on the second row 
the ON-model with a — 0.1, on the third row the PH-model with ks = 0.2. 


8 Conclusions 

In this paper we derived nonlinear partial differential equations of first or¬ 
der, i.e. Hamilton-Jacobi equations, associated to the non-Lambertian re¬ 
flectance models ON-model and PH-model. We have obtained the model 
equations for all the possible cases, coupling vertical or oblique light source 
with vertical or oblique position of the observer. This exhaustive description 
has shown that these models lead to stationary Hamilton-Jacobi equations 
with a same structure and this allows for a unified mathematical formula¬ 
tion in terms of fixed point problem. This general formulation is interesting 
because we can switch on and off the different terms related to ambient, 
diffuse and specular reflection in a very simple way. As a result, this gen¬ 
eral model is very flexible to treat the various situations with vertical and 
oblique light sources. Unfortunately, we have observed that none of these 
models is able to overcome the typical concave/convex ambiguity known 
for the classical Lambertian model. Despite this limitations, the approach 
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Figure 21: Experiments on real images: input images and the recovered 
shapes obtained by our approach. On the first row: urn reconstructed by 
the L-model with uj — (0, 0,1). On the second row: rabbit reconstructed by 
the ON-model with a = 0.2, oj — (0,0,1), V = (0,0,1). On the third row: 
Beethoven reconstructed by the PH-model with uj — (0.0168,1.198, 0.9801), 
ks = 0.2. 

presented in this paper is able to improve the Lambertian model that is 
not suitable to deal with realistic images coming from medical or security 
application. The numerical methods presented here can be applied to solve 
the equations corresponding to the ON-model and the PH-model in all the 
situations although they suffer for the presence of several parameters re¬ 
lated to the roughness or to the specular effect. Nonetheless, looking at the 
comparisons with other methods or models shown in this paper, one can 
see that the Semi-Lagrangian approach is competitive with respect to other 
techniques used, both in terms of CPU time and accuracy. As we have seen. 
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in the complex nonlinear PDEs associated to non-Lambertian models the 
parameters play a crucial role to obtain accurate results. In fact, varying 
the value of the parameters it is possible to improve the approximation with 
respect to the classical L-model. We can also say that for real images the 
PH-model seems easier to tune, perhaps because we need to manage less 
parameters. 

Focusing the attention on the tests performed with an oblique light source, 
we have to do some comments that are common to the PH-model and the 
ON-model. Several terms appear in these models and each of them gives 
a contribution to the roundoff error. Note that the accumulation of these 
roundoff errors makes difficult in the oblique case to obtain a great accu¬ 
racy. A possible improvement could be the use of second order schemes, 
that release the link between the space and the time steps which charac¬ 
terizes and limits the accuracy for first order schemes. Another interesting 
direction would be to extend this formulation to other reflectance models 
(like e.g. the Ward’s model) and/or to consider perspective projection also 
including an attenuation term which can help to resolve the concave/convex 
ambiguity or considering more than one input image with non-Lambertian 
models, that solves the well-known ambiguity (see 131 for a first step in this 
last direction considering the Blinn-Phong model). These directions will be 
explored in future works. 
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